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Abstract 

Polarization in ferroelectrics, described by the Landau- Ginzburg Hamiltonian, is 
considered, based on a multi-dimensional Fokker-Planck equation. This formulation 
describes the time evolution of the probability distribution function over the polar- 
ization field configurations in the presence of a time dependent external field. The 
Fokker-Planck equation in a Fourier representation is obtained, which can then be 
solved numerically for a finite number of modes. Calculation results are presented for 
one and three modes. These results show the hysteresis of the mean polarization as 
well as that of the mean squared gradient of the polarization. 

Keywords: Ferroelectric, Landau-Ginzburg Hamiltonian, Fokker-Planck equation, 
polarization hysteresis. 
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1 Introduction 

The stochastic description of collective phenomena such as phase transitions is a key 
theme in solid-state physics. Problems of this kind are nontrivial and are usually solved 
by means of perturbation theory [USE]. Here the kinetics of polarization switching in 
ferroelectrics is studied, taking into account the spatio-temporal fluctuations of the polar- 
ization field given by the Langevin and multi-dimensional Fokker-Planck equations [U [5] . 
Simplified approaches are well known U ||, for which the static distribution of the 
polarization is found by minimizing the Landau-Ginzburg type free energy functional. 
The dynamics and fluctuations of the polarization are usually introduced via equations 
of mean-field type representing either an approximation or a solution of the appropriate 
mean- field model [9] . The multi-dimensional Fokker-Planck equation discussed here allows 
us to consider unlimitedly many nontrivial degrees of freedom for the correlated collective 
fluctuations of the polarization field, which are lost or roughly treated as uncorrelated in 
the mean-field approximation. This problem has been studied previously [101 111] by means 
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of the perturbative Feynman diagram technique. Such an approach, however, only allows 
the description of the ferroelectric response with respect to an infinitely small external 
field, rather than polarization reversal (switching) and hysteresis. The latter requires a 
non-perturbative treatment |12[ [T3| [Tl] which, however, is usually done by including some 
kind of mean-field approximation. In [15 1 116 1 [T7] the polarization reversal has been studied 
based on the classical nucleation theory. At the initial stage, the formation of domains of 
new phase is considered as one-step Markov process, in which case the domain size can 
increase or decrease by one unit at a time moment. The obtained in this way distribution 
over the sizes of small fluctuating domains (or nuclei) is then used to describe the later 
stage of the domain growth and the polarization reversal. Despite of some simplifications 
(neglecting nontrivial spatial correlations as well as complicated merging and splitting 
processes) it gives realistic results [151 [T7] . An explicit consideration of the nucleation and 
growth of domains lies also in the basis of more recent simulations of ferroelectricity in 
polycrystals and single crystals |T8l H9] , as well as of combined ferro- and piezo-electric 
effects [201 [21]. These ideas are applied also to describe experimentally observed features 
in specific materials [22, 23J. In the present work the Fokker-Planck equation in a Fourier 
representation, which is suitable for a non-perturbative numerical approach avoiding ap- 
proximations of the mean-field type, is derived. In principle, it allows to include many 
nontrivial degrees of freedom, lost in the known simplified approaches, as regards the spa- 
tial correlations and fluctuations of the polarization. In practical calculations of this kind, 
however, only some fluctuation degrees of freedom, that is only a few Fourier modes, can 
be implemented due to computation time and memory constraints. Nevertheless, such 
limited calculations can reproduce some important qualitative features and provide test 
examples to verify possible approximations which would allow the treatment of realistic 
numbers of modes. 

2 Basic equations 

We consider a ferroelectric with the Landau-Ginzburg Hamiltonian 

H = J (|p 2 (x) + ^P 4 (x) + |(VP(x)) 2 -A(x,t)P(x,t)) ( ix, (1) 

where P(x, t) is the local polarization and A(x, t) is the time-dependent external field. 
The only configurations of the polarization which are allowed are those corresponding to 
the cut-off k < A in Fourier space with A = n/a, where a is the lattice constant. The 
Hamiltonian ([I]) can be approximated by a sum over discrete cells, where the size of one 
cell can be larger than the lattice constant. Such cells are small domains with almost 
constant polarization. Thus the Hamiltonian of a system with volume V is 

H = AV E (f p2 w + i p4 w + |( vp ( x )) 2 - A ( x < W x ' *)) ■ ( 2 ) 

where AV = V/N is the volume of one cell and the coordinates of the centers of the cells 
are given by the set of discrete d-dimensional vectors x G R d . The stochastic dynamics of 
the system is described by the Langevin equation 

dH 
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where £(x, t) is white noise, i.e. 

(e(x,i)e(x / ,t , )) = 27^x,x'^-t')- 

In the case of Gaussian white noise, the probability distribution function 

/(P(x 1 ),P(x 2 ),---,P(x JV ),t) 
is given by the Fokker-Planck equation [5] 



(4) 
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At equilibrium the flux vanishes, which corresponds to the Boltzmann distribution / oc 
exp(-H/9) with 9 = k B T. 

Assuming periodic boundary conditions, we consider the discrete Fourier transform 

(7) 



p( X ) = N-^j^Pkj 

k 

Pk = N-^^P^e 



ikx 



The Fourier amplitudes are the complex numbers P k = P k + iP k . Since P(x) is real, 
P'_ k = P k and P" k = — P k hold. It is supposed that the total number of modes N is an odd 
number. This means that there is a mode with k = and modes with ±ki, ±k2, . . . , ±k m , 
where m = (N — l)/2 is the number of independent non-zero modes. 
The Fokker-Planck equation for the probability distribution function 



/ — / (-Po, -Pfci > P k 2 ' • • • ' P k rn i Pk! ) Pk 2 ' • • • ' p L m > *) 
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is then 
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where Pg = Po, ^ is the set of m independent non-zero wave vectors and includes k = 0. 
Here the Fourier transformed Hamiltonian is given by 



H 



Ay 



f^E( a + ck2 ) i p ki 2 +^ 1 E p^p^a 

\ k ki+k 2 +k 3 +k 4 =0 

-^A_ k (i)P k ) . 



(10) 



Some of the variables in Eq. (fTUj) are dependent according to P'_ k = P k and P" k = 
— P k - Here A k (t) = A k (t) + »A k (t) is the Fourier transform of A(x, t). The Fokker-Planck 
equation Q can be written explicitly as 

ldf 

7 at 



= E apT { AV / [(a + ck 2 ) P k + /?5 k - A k (t)] + | (1 + 5 k , ) JO 
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where 



= N- 1 {Wka-<«}. (12) 

ki+k2+k3=k 

S£ = £ + (13) 

k!+k2+k 3 =k 



3 Spatially homogeneous case 

The simplest case is a spatially homogeneous polarization for which only the k = mode 
is retained in (jlip with a spatially homogeneous external field \(x,t) = Xo(t) = A$m(cjt). 
In this case we have 

^%=M vflaPo+f)FS - AM " t)]+ ''£k} ■ (i4) 

The mean value of the zero-th Fourier amplitude 

oo 

Po= J Pof(Po,t)dP (15) 

— oo 

is the mean polarization P in this case, as follows from ©. 



4 Quasi one-dimensional case with homogeneous 
external field 

We shall further consider a quasi one-dimensional case, for which a three-dimensional 
ferroelectric sample is stretched out in the x-direction, i.e. L x 3> L y and L x 3> L z hold 
for the linear sizes. In this case we assume that the polarization, as well as the external 
field, depend only on the x coordinate. This means that the wave vectors also have 
only one non-vanishing component, which is a scalar quantity k n = (2tt/L x ) ■ n, where 
n = 0, ±1,±2, ... , ±m. From now on, we shall omit the vector notation in this quasi 
one-dimensional case. 

As a first step, we shall include only one (m = 1) independent non-zero wave vector 
ki = 2tt/L x (in total N = 3 wave vectors k = —kL,0,ki) and a homogeneous external 
field 

\(x,t) = ^=X (t) = ^=Asm(ut). (16) 

Furthermore, we shall assume that the probability distribution function in real space is 
translation invariant at the initial time. Due to the translational symmetry of the model, 
translation invariance then also holds at all later times. In the Fourier representation, this 
means that the probability distribution function depends on the modulus of , but not 
on its phase. Thus we have 

/ = / (P , P' kl , Pi , t) = nP ^\ P p^\ t) ' (17) 
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where f (Po, \ P kl \,t) is the probability density in the (Pq, \ P kl |) space and which obeys 
the Fokker-Planck equation 
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Since / is finite, / vanishes at | P kl | = 0. The physical boundary conditions correspond 



to zero flux at the boundaries Pn = ±oo, 



and Pi 



ki 



oo. An appropriate 



initial condition has to be chosen which fulfils these relations, e.g. / (i-b> I Pki |> 0) oc| P kl 
exp (-aoP 2 - ai | P kl | 2 ). 



5 Numerical results 



The Fokker-Planck Eq. (lllj) has been solved numerically in the quasi-one dimensional case 
for which the wave vector has only one component k = k x , in the spatially homogeneous 
approximation with only the n = (k x = 0) mode (in total N = 1 modes), as well as in 
the next higher order approximation for which the spatial distribution of the polarization 
field is taken into account by including n = 0, =Fl (in total N = 3 modes) with k n = k^ n. 
At first, rewriting the Eq. in the following form: 



ldf d i d i d a 

- 1 m=dpr o J ^dPi i J ^dPj i J ^ (19) 



where (for a < 0): 



r n = AV [(M - ck n 2 ) PI - (3SI + XtJt)) f - rZ 9 --^, (20) 

fori 

< = { l + 5kn 'ZZ, l , - = 0,L (21) 
Then, by using a special exponential-type substitution for the distribution function: 

/ = Wt n exp jjf *" ^ [(M - ck n 2 ) p - (3St n + XIM dpj , (22) 

where Wt is a normalization function and po is a real number (similarly as in [23], it can 
be shown that W k and po should not affect final coefficients of a difference scheme), a 
monotone, exponential difference scheme has been developed: 

(A fo) = £ ±- (Af)' + ±. Wj d n - (Qi)' H +1 + 

n=0 i' n 4 
n^lK'i ^ ^ " 7 Tl+l 

< i* < 
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with coefficients: 



(exp f (7/),_ 1/2e& 1 -l j 



exp 
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= ^ ((-I«l + ^„ 2 ) ^ + PSt n - X kn (t l+1 )) h it , (27) 

where i = (io>*ij*i) an d is unit vector in the direction, I is the time index, h is the 
polarization space step and r is the time step. 

Elaborated difference scheme (123p - (127p has a first order truncation error in time and 
a second order truncation error in space, when \rj\ — > 0. In the limit |ry| — > oo, so that 
transport is advection dominated, when 6 — > 0, the difference scheme (j23|) - (|27p has first 
order truncation errors in time and space. 

For one dimensional case, when N = 1, difference scheme (|23|) - (|27p coincides with 
elaborated in [25]. As it was shown [25], for \rj\ — > and \rj\ — > oo one dimensional 
exponential scheme becomes absolutely stable. Here in order to control the precision 
(stability) of the calculations we adopted a linear prognosis criterion. Accordingly, we 
required that the linear prognosis at each time step should not deviate from the numerical 
solution to a given precision, which varied from si = 2 ■ 10 -4 to 6f = 5 • 1CP 4 . If this 
precision criterion was met, then the time step was kept constant. The time step was 
increased if the linear prognosis deviation was less than ei, and it was decreased if the 
linear prognosis deviation was greater than ej. 

Elaborated absolute monotone scheme (|23p - (|27p for the Fokker-Planck Eq. (jlip has 
been implemented on parallel processing for one and three polarization space dimensions, 
when N = 1 and N = 3. A parallel code version of program is written by implementing 
MPI programming technology in FORTRAN and using ScaLAPACK (SUN S3L) linear 
algebra solver. 

In contrast to the quasi-one dimensional case discussed in Sec. HI here we consider 
homogeneous as well as non-homogeneous external fields based on (jlip as a generic 
equation. This equation reduces to (|14p for N = 1. In the case for which the field 
is spatially homogeneous it is given by Ao(t) = Asin(ujt) = A'yfN sin(a;t), the other 
Fourier components being zero, whereas the non-homogeneous field is chosen such that 
Ao(i) = A' fci (i) = Aj£ At) = A't/N sin(ut). In the coordinate representation, this corre- 
sponds to X(x,t) = ^4'sin(u;i) in the homogeneous case and to 

X(x, t) = A' sin(cjt) 1 + 2V2 cos(A;ix + vr/4) 

in the non-homogeneous case. Note that, according to ©, the mean polarization is given 
by P = Po/y/N. We shall also compute the mean-squared gradient of the polarization 
(VP) 2 , which is given by 

J^ = N- 1 Y J ^ 2 \P^ = N- 1 Y J k 2 (P^ + Pf) (28) 
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in general. In our quasi one-dimensional case the sum over k reduces to the sum over 
k = -k L ,0,k L . 

In the Eq. (|lip expressions for S were given by Sq = Pq for N = 1, and for N = 3 
these are: 

So = i(Po) 3 + 2Po((^ 1 ) 2 +(^' 1 ) 2 ), (29) 
S' kl = (pi)\pi((Po)\(Pk 1 ) 2 ), (30) 



S k \ = {P kl ) + P k \{{ P o) + [P kl ) ) ■ (31) 

The numerical calculations have been performed within the interval [—3.5, 3.5] for each 
Fourier amplitude (Po f° r N = 1 and Pq, P' k , P{! for N = 3) with the boundary condition 
/ = at the borders of this domain. The parameter values are 7 = V = /3 = 1, L x = 7, 
a = -1,9 = 0.05, A' = 0.5 (i.e. A = 0.5-/2V) and cj = 10~ 3 . The initial condition 



/(P ,0) = ^^exp 



\ 



2a 2 



(32) 



has been used in the spatially homogeneous (JV = 1) case and has been modified to 

(33) 



/(^P^P^OH^^exp 



( ( p*-p) 2 + p' kl 2 + p'C 



V 



2a 2 



in the other case of N = 3 modes, with a = 0.3, P = — 1 for Po < and P = 1 for Po > 0. 

The calculated mean polarization in the coordinate space P = Po/yN, which depends 
on the homogeneous external field represented by \o(t), forms a hysteresis loop, as shown 
in Fig. [H where the results for N = 1 (dot-dashed curve) and N = 3 (solid curve) are 
compared. In the second case the spatial distribution of the polarization, generally, is 
nonhomogeneous. The insertions in Fig. [T] illustrate qualitatively some of possible or 
typical polarization distributions at different locations on the hysteresis loop. As can 
be seen, the inclusion of iV = 3 modes only slightly changes the behaviour of the mean 
polarization as compared with the homogeneous approximation with only the zero mode 
k = included (N = 1 case). The hysteresis loop becomes narrower for N = 3. This 
is clear physically, since N = 3 additional degrees of freedom for the fluctuations in the 
polarization make reversal of the polarization easier. This phenomenon is related to the 
known Landauer paradox. Namely, the observed coercive field in real ferroelectrics is 
much smaller than that calculated assuming a spatially homogeneous switching of the 
polarization. It is well known that switching in ferroelectrics is a process driven by the 
nucleation and growth of domains of different polarization [15j [T71 QjJl [22]. Hence, the 
above mentioned effect of narrowing the hysteresis loop should be even much stronger 
for appropriate system parameters and larger number of Fourier modes N, allowing to 
describe the domain structure realistically. 

The effect of non-homogeneity of the external field on the polarization hysteresis can 
be seen in Fig. [21 where numerical results for the homogeneous (solid curve) and the 
above discussed non-homogeneous (dashed curve) fields are compared in the three mode 
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-0.6 XJt)/A 0.6 



Figure 1: The polarization hysteresis: the mean polarization P vs. normalized external 
field \o(t)/A calculated numerically for the parameter values j = V = /3 = l, a = — 1, 
6 = 0.05, A = and u = 10 -3 . The solid and dot-dashed lines represent the 

numerical results for N = 3 and N = 1 Fourier modes, respectively. The insertions 
illustrate possible spatial distributions of the polarization at different locations on the 
hysteresis loop. 




Figure 2: The polarization hysteresis: the mean polarization P vs. \q{€)/A calculated 
numerically for N = 3 Fourier modes. The solid and the dashed lines represent the results 
for a homogeneous and non-homogeneous external field, respectively. 



(VP) 2 
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0.02 
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0.01 
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Figure 3: The mean squared gradient of the polarization (VP) 2 vs. Xo(t)/A calculated 
numerically for N = 3 Fourier modes for the case of a homogeneous external field. The 
pieces of the hysteresis loop are numbered in order of increasing time. 

(N = 3) approximation. The maximum values of | P | are somewhat smaller for the 
non-homogeneous field. This is related to the fact that the non-homogeneous field 



changes sign depending on the coordinate x in such a way that it is positive in 61.5% of 
the volume and negative in the rest, or vice versa. This means that the polarization, with 
a large probability, has an opposite sign in the smaller part of the volume, which reduces 
the modulus of the mean polarization. 

In contrast to the spatially homogeneous approximation with N = 1, inclusion of the 
non-zero modes allows us to evaluate not only the mean polarization, but also the mean 
squared gradient of the polarization (j28H . which is a measure of the spatial inhomogeneity 
of P(x, t). The hysteresis loop for this quantity, calculated for N = 3 for the homogeneous 

external field, is shown in Fig. An interesting feature is that (VP) 2 has maxima at times 
which roughly correspond to those at which the fastest switching of the mean polarization 
takes place. This provides evidence that the polarization switching quite often, with a 
remarkable probability, is accompanied by formation of a spatially inhomogeneous struc- 
ture. The insertions in Fig. [1] illustrate such a scenario. By including more Fourier modes, 
this eventually could be identified with a domain-like structure. This is the most probable 
scenario expected from the known studies of the polarization reversal as a nucleation and 
domain growth process (see, e. g., [151 [13 HH] and references therein). 

In Fig. HI the (VP) 2 hysteresis loop is plotted for the non- homogeneous external field. 
The maxima related to the polarization switching are observed in this case too. Contrary to 
the case for the homogeneous field, which tends to order the polarization in one direction, 
the sign alternating non-homogeneous field forces the formation of a non-uniform spatial 
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Figure 4: The mean squared gradient of polarization (VP) 2 vs \o(t)/A calculated numer- 
ically for N = 3 Fourier modes for the case of the non-homogeneous external field. Pieces 
of the hysteresis loop are numbered in order of increasing time. 

distribution in the polarization P(x,t). It can be seen from Figs. [3] and H] that the mean 
squared gradient of P(x, t) decreases at the largest absolute values of the field if the field 
is homogeneous (Fig. [3]) and, remarkably, increases for the non-homogeneous field (Fig. 2]). 

6 Conclusions 

In the present work a multi-dimensional Fokker-Planck equation has been derived which 
describes the kinetics of polarization switching in a ferroelectric in the presence of an ex- 
ternal field. The probability distribution function for this equation depends on a set of 
Fourier amplitudes. Calculations were performed for a spatially homogeneous approxima- 
tion retaining only the zero mode k = 0, as well as for three Fourier modes {k = —hz, 0, k^) 
for both homogeneous and non-homogeneous external fields. The hysteresis of the mean 
polarization P and of the mean squared gradient of the polarization (VP) 2 were calculated 
and compared. In particular, it was found that the P hysteresis loop becomes slightly nar- 
rower when more Fourier modes are included. This loop is qualitatively the same as that 
observed in real ferroelectrics. The mean squared gradient of polarization is a measure 
of its inhomogeneity. The hysteresis loop of (VP) 2 provides evidence that polarization 
switching is accompanied by an increased spatial inhomogeneity of P(x,t), as expected 
from the known theoretical approaches describing the polarization reversal as a nucleation 
and growth of domains. 

Although only a few Fourier modes were included in the actual calculations, the prob- 
lem was treated non-perturbatively and without the mean-field approximation. More 
precisely, our calculations were based on exact equations for the probability distribution 
function in the case for which the Hamiltonian (llOj) contains a given number of modes. 
The present results thus can be used to test various possible approximations. 
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